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We  develop  and  evaluate  a  high-order  discontinuous  Galerkin  method 
for  the  solution  of  the  shallow  water  equations  on  the  sphere.  To  overcome 
well  known  problems  with  polar  singularities,  we  consider  the  shallow  water 
equations  in  Cartesian  coordinates,  augmented  with  a  Lagrange  multiplier 
to  ensure  that  fluid  particles  are  constrained  to  the  spherical  surface. 

The  global  solutions  are  represented  by  a  collection  of  curvilinear  quadri¬ 
laterals  from  an  icosahedral  grid.  On  each  of  these  elements  the  local 
solutions  are  assumed  to  be  well  approximated  by  a  high-order  nodal  La¬ 
grange  polynomial,  constructed  from  a  tensor-product  of  the  Legendre- 
Gauss-Lobatto  points  which  also  supplies  a  high-order  quadrature.  The 
shallow  water  equations  are  satisfied  in  a  local  discontinuous  element  fash¬ 
ion  with  solution  continuity  being  enforced  weakly. 

The  numerical  experiments,  involving  a  comparison  of  weak  and  strong 
conservation  forms  as  well  as  the  impact  of  over-integration,  confirm  the  ex¬ 
pected  high-order  accuracy  and  the  potential  for  using  such  highly  parallel 
formulations  in  numerical  weather  prediction. 


Key  Words:  discontinuous  Galerkin  method,  liigh-order,  icosahedral  grid,  shallow  water 
equations,  spectral  element  method,  spherical  geometry. 


1.  INTRODUCTION 

The  majority  of  current  climate  and  numerical  weather  prediction  (NWP)  mod¬ 
els,  e.g.,  the  operational  NWP  models  developed  by  the  European  Center  for 
Medium-Range  Weather  Forecasting  (ECMWF),  the  National  Center  for  Environ¬ 
mental  Prediction  (NCEP),  and  the  Naval  Research  Laboratory  (NRL),  are  based 
on  globally  defined  spectral  methods  [11,  14,  24,  25].  While  these  methods  continue 
to  be  essential  tools  for  NWP,  their  limitations  are  beginning  to  emerge.  For  one, 
the  fixed  global  grid  makes  adaptive  solution  techniques  very  complex  if  possible  at 
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all.  Furthermore,  the  recent  paradigm  shift  in  large-scale  computing  from  vector 
to  distributed-memory  computing  platforms  has  exposed  problems  with  achieving 
efficiency  due  to  the  global,  inter-processor,  all-to-all  communication  needed  in  the 
spectral  transform. 

Such  considerations  have  stimulated  research  into  parallel  methods  and  poten¬ 
tially  adaptive  methods,  suitable  for  NWP,  yet  maintaining  the  accuracy  of  the 
global  spectral  approach  while  overcoming  its  inherent  limitations. 

Continuous  spectral  element  methods  have  recently  been  proposed  for  future 
climate  [26]  and  NWP  [9]  models.  In  these  methods  the  solution  is  first  constructed 
in  a  local  element-by-element  manner  and  then  a  portion  of  this  local  solution  is 
distributed  to  the  global  grid  points  shared  by  adjacent  elements.  This  global 
assembly  from  the  elements  to  the  global  grid  points  is  what  enforces  continuity 
and  it  is  this  requirement  which  reduces  the  locality  of  the  method.  Although 
continuous  spectral  element  methods  parallelize  quite  well  [27]  their  insistence  on 
continuity  makes  them  cumbersome  to  employ  in  either  a  non-conforming  approach 
[17]  or  in  an  adaptive  solution  strategy. 

In  this  paper  we  introduce  a  nodal  high-order  discontinuous  Galerkin  method  for 
geophysical  flows  on  the  sphere.  Like  continuous  spectral  element  methods,  discon¬ 
tinuous  Galerkin  methods  (DGM)  can  be  constructed  to  have  high-orcler  accuracy, 
while  maintaining  a  large  degree  of  locality,  hence  enabling  high  parallel  perfor¬ 
mance  and  adaptive  solution  procedures.  The  locality  of  these  methods  ensures 
that  they  can  be  used  with  any  type  of  grid,  e.g.,  unstructured  and  non-conforming 
if  needed.  In  this  paper  we  shall  demonstrate  this  by  using  unstructured  icosahedral 
grids. 

As  a  first  step  towards  the  construction  of  a  fully  3D  atmospheric  model  we 
demonstrate  the  efficiency  and  accuracy  of  the  nodal  high-order  DGM  by  solving 
the  shallow  water  equations  on  the  sphere.  The  shallow  water  equations  contain  all 
of  the  horizontal  operators  required  in  an  atmospheric  model  and  thus  represent  a 
good  first  test  for  newly  proposed  methods  for  atmospheric  models. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Sec.  2  we  introduce  the 
spherical  shallow  water  equations  and  discuss  the  reasons  for  using  a  Cartesian 
grid.  This  sets  the  stage  for  Sec.  3  which  introduces  the  numerical  scheme  and 
discusses  in  detail  the  spatial  curvilinear  representation  of  the  solution  as  well  as 
the  formulation  enabling  one  to  satisfy  the  equations  in  a  discontinuous  fashion.  We 
also  discuss  the  temporal  time-stepping  scheme  and  the  approach  taken  to  generate 
a  suitable  grid  on  the  sphere.  In  Sec.  4  we  demonstrate  the  accuracy,  efficiency,  and 
robustness  of  the  complete  scheme  for  the  solution  of  benchmark  problems  for  the 
spherical  shallow  water  equations.  Section  5  contains  a  few  remarks  and  outlines 
natural  extensions  of  the  work  presented  here. 

2.  SPHERICAL  SHALLOW  WATER  EQUATIONS 

The  spherical  shallow  water  equations  in  conservation  form  are  given  as 


where 


-^  +  \/-F(q)  =  S(x,q) 


(1) 
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q=[tp,<pu,tpv,tpw]T  , 

represents  the  state  vector,  q ,  composed  of  the  geopotential  height,  p,  and  the  three 
Cartesian  velocity  components,  ( u,v,w ),  all  being  a  function  of  x  €  R'5  and  time, 
t.  The  units  for  p  and  u  are  (m/s)2  and  m/s,  respectively. 

The  flux,  F(q),  takes  the  form 


F(q) 


ipu 

ipv 

(fW 

Lpu2  +  \y2 

i  + 

( puv 

2  i  1  2 

3  + 

(fUW 

( puv 

< pvz  + 

( pvw 

(fUW 

(fVW 

_  pw 2  +  \p 2  . 

(2) 


where  ( i,j,k )  represent  the  Cartesian  unit  vectors.  The  source  term,  S(x,q)7  in 
Eq.(l),  acting  only  on  the  momentum  equations,  is  given  as 

.  2f lip 

S(x7  q)  =  -——x  xu-  p\7ps  +  px  . 

Rz 

Here  the  first  term  accounts  for  the  Coriolis  force,  with  R  —  6.371  x  106  m  being 
the  radius  and  O  =  7.292  x  10-5  rad/s  the  angular  frequency  of  the  earth,  while 
the  second  contribution  models  the  effects  of  a  variable  surface  height  through  the 
surface  potential,  (ps.  The  last  term  is  a  Lagrange  multiplier,  the  specification  of 
which  we  shall  return  to  shortly. 

Equation  (1)  is  derived  from  the  incompressible  and  inviscid  Navier-Stokes  equa¬ 
tions  by  vertically  integrating  the  mass  to  yield  the  geopotential  height  equation 
(for  further  details  see  [22]).  Contrary  to  most  other  work  on  the  numerical  solu¬ 
tion  of  Eq.(l)  on  a  spherical  shell,  we  shall  not  recast  it  in  spherical  coordinates 
but  rather  maintain  the  Cartesian  coordinates.  The  main  motivation  for  doing  so, 
albeit  at  the  expense  of  introducing  an  additional  momentum  equation,  is  to  avoid 
the  problems  associated  with  the  polar  singularity.  For  a  spherical  shell,  described 
by  the  coordinates  (A,  9),  of  radius  R  the  divergence  of  a  vector  field,  F  =  f  X  +  g9, 
is  given  as 


V  •  F  = 


1 

R  cos  9 


'd£ 

d\ 


dg  cos  9 
89 


At  the  poles,  i.e.,  9  =  ±7r/2,  this  is  a  source  of  numerical  problems,  caused  by  the 
specific  formulation  rather  than  the  nature  of  the  shallow  water  equations  and  its 
solutions.  While  the  use  of  a  local  Cartesian  coordinate  system  has  been  used  to 
overcome  these  problems  in  the  past  [26]  we  have,  guided  by  the  results  of  previous 
work  [9],  chosen  to  maintain  the  Cartesian  formulation  everywhere. 

To  ensure  that  the  fluid  particles  remain  on  the  spherical  shell,  we  require  that 
the  fluid  velocity  remains  perpendicular  to  the  position  vector,  i.e.,  x  ■  u  =  0,  which 
yields  the  Lagrange  multiplier 


1 

If 


px  ■  Vps  +  X  ■ 


p(x,q) 


(3) 
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Here  F  represents  the  parts  of  the  flux,  Eq.  (2),  associated  with  the  momentum 
equations,  i.e.,  the  last  three  fluxes.  Note  that  the  term  associated  with  the  Coriolis 
term  vanishes  identically  at  this  exact  level  while  this  may  not  be  case  at  the 
discrete  level.  As  no  constraint  is  needed  on  the  geophysical  potential,  essentially 
representing  the  local  mass,  the  multiplier  is  needed  in  the  momentum  equations 
only. 


3.  THE  NUMERICAL  SCHEME 

In  developing  the  numerical  scheme  for  the  solution  of  Eq.  (1)  we  shall  split  the 
discussion  into  a  treatment  of  the  spatial  discretization  and  the  approximation  of 
the  resulting  semi-discrete  approximation. 


3.1.  The  Spatial  Approximation 

The  discussion  of  the  spatial  approximation  scheme  involves,  as  does  any  formula¬ 
tion  of  a  scheme  for  solving  partial  differential  equations,  attention  to  the  questions 
of  how  to  represent  the  solution  as  well  as  in  which  way  the  equations  are  required 
to  be  satisfied.  In  the  following  we  address  these  two  issues  in  more  detail. 


3.1.1.  Representing  the  Solution  and  Basic  Operations 

Initially,  we  assume  that  the  computational  domain,  S,  i.e.,  a  spherical  shell,  is 
covered  by  K  non-overlapping  curvilinear  quadrilaterals,  D,  such  that 

S  =  (J  Dfe  . 

k= 1 


The  construction  of  this  sphere  covering  is  not  entirely  trivial  and  we  shall  return 
to  this  problem  in  Sec.  4.  For  now,  however,  we  simply  assume  its  existence. 

To  enable  the  efficient  and  accurate  computation  of  operations  such  as  differ¬ 
entiation  and  integration,  we  introduce  a  nonsingular  mapping,  x  =  4'(£),  which 
connects  the  local  physical  coordinates,  x  =  ( x,y,z ),  defined  on  D  with  a  refer¬ 
ence  system  £  =  (£,77,  £),  defined  on  the  local  element  such  that  (£,77)  lies  on  the 
spherical  surface.  Thus,  £  represents  the  spherical  radius  vector  itself,  i.e.,  £  con¬ 
stant  corresponds  to  a  shell  of  a  constant  radius.  For  simplicity  we  assume  that 
(£,  77)  G  [—1,  l]2  on  each  element,  see  Fig.  1. 

Associated  with  the  local  mapping,  4',  is  the  transformation  Jacobian,  J  = 
and  the  determinant 


\J  I 


dx  _  _  dx  dx 

8Tg  ’  G=a?  xat 


where  G  represents  the  surface  conforming  component  of  the  mapping  (see  [9]  for 
further  details). 

The  mapping  also  supplies  the  local  metric,  V£  and  V77,  as  well  as  the  local 
normal  vectors  at  the  surface  of  the  element.  Indeed,  as  illustrated  in  Fig.  1,  the 
normals  are  given  as 


fi  =  rj 


Wv\ 


,  n  =  £ 


V£ 


|V£| 


77= ±1 


£=±1 
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FIG.  1.  The  geometry  of  the  mapping  and  the  associated  metric  information. 

for  the  sides  1/3  and  2/4,  respectively. 

With  this  in  place  we  can  now  focus  on  the  local  element-wise  representation  of 
the  solution,  q,  and  the  approximation  of  operations  such  as  differentiation  and 
integration.  For  simplicity,  we  assume  £  to  be  unity  in  what  remains  and  denote 
£  =  (£,  q)  unless  clarification  is  deemed  necessary. 

The  simple  structure  of  the  standard  element,  I,  spanned  by  £  £  [—1,  l]2,  makes 
it  natural  to  represent  the  local  solution  by  an  TVth  order  polynomial  in  £  as 

(N+l)2 

Qn(x)  =  ClN(xk)Lk{x)  , 

fc= 1 

where  xj.  represents  (./V+l)2  grid  points  and  Lk(x)  reflects  the  associated  multivari¬ 
ate  Lagrange  interpolation  polynomial.  The  logical  square  structure  of  I  simplifies 
matters  further  in  that  we  can  express  the  Lagrange  polynomial  by  a  tensor-product 
as 


Lk(x)  =  ,  (4) 

where  i,  j  =  0, ...,  N. 

While  many  choices  of  the  grid  points,  (£,,  rjj ),  are  possible,  it  is  natural  to  choose 
the  Legendre-Gauss-Lobatto  points,  given  as  the  tensor-product  of  the  roots  of 


(1-£2)PM£)=0 


where  TW(£)  is  the  ATtli  order  Legendre  polynomial.  This  choice  is  natural  as  these 
points  are  endowed  with  a  Gaussian  quadrature  rule  which  shall  become  useful 
shortly.  With  this  choice,  the  one-dimensional  Lagrange  polynomials,  /ii(£),  also 
known  as  the  Legendre  cardinal  functions,  become 


hi(  0  = 


i 

N(N  +  1)  (£-6)^(6) 
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and  likewise  for  hj  ( ? )) . 

The  choice  of  the  Legendre-Gauss-Lobatto  points  enables  the  straightforward 
approximation  of  element-wise  integrals,  i.e. , 


[  q(x)dx  =  [ q{£)J(Q 

Id  J  i 


d( 


N 

i,j= 0 


where  J  represents  the  local  Jacobian  for  the  transformation  between  D  and  I,  and 
u>l  and  wj  are  the  Gaussian  quadrature  weights, 


2 

N{N  +  1)  \Pn{^)) 


associated  with  the  one  dimensional  Legendre-Gauss-Lobatto  quadrature.  We  recall 
that  if  q  J  is  a  polynomial  of  at  most  degree  2N  —  1  in  each  local  coordinate,  the 
quadrature  is  exact. 

The  evaluation  of  surface  integrals,  made  particularly  simple  by  the  natural  sep¬ 
aration  between  interior  and  edge  (side)  nodes  in  the  nodal  formulation  considered 
here,  follows  the  same  line  of  thinking,  i.e., 


®  q(x)  dx 
JSD 


l  mmdt 

JS I 
N 

X  [q{Zi,  -1  -1)  +  «(&,  1  )j&,  1)]  4 

i= 0 

N 

j= o 


3.1.2.  Satisfying  the  Equation 

With  the  local  representation  of  the  solutions  and  the  central  operations  in  place, 
we  shall  now  proceed  to  consider  the  question  of  how  to  satisfy  the  equation.  We 
assume  that  the  solution,  q,  to  Eq.(l)  is  represented  locally  by  high-order  polyno¬ 
mials,  qN,  defined  on  the  curvilinear  quadrilateral,  D,  and  require  that  the  equation 
be  satisfied  element-wise  in  the  following  discontinuous  Galerkin  way 


+  V  •  Fn 


Lk(x)  dx 


Lk(x)h  ■  [F N  -  F*N]  dx  ,  (5) 


where  Lk(x)  (k  G  [1, ...,  (N  +  l)2])  is  the  local  polynomial  basis,  Eq.(4).  We  have 
also  introduced  the  polynomial  representation  of  the  flux,  Fjy  =  Fpf(qN),  and  the 
source,  Sn  =  S>N(x,qN),  as 


(AT+l)2  (iV+l)2 

FN{qN)=  X  F{ciN(xk))Lk(x)  ,  SN(x,qN)=  X  S{xklqN{xk))Lk(x)  . 

k=  1  k= 1 
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The  numerical  flux,  F*N,  in  Eq.(5)  shall  be  discussed  in  detail  shortly.  Prior  to  that, 
however,  a  few  remarks  concerning  Eq.(5),  are  in  order.  First  of  all  we  note  that 
the  interface  conditions,  introduced  into  the  formulation  through  the  numerical 
flux,  are  enforced  only  weakly,  i.e.,  the  solution  is  in  general  discontinuous.  As 
we  shall  see,  however,  this  does  not  impact  the  accuracy  as  the  size  of  the  jump 
vanishes  to  the  order  of  the  interior  approximation.  Furthermore,  the  discontinuous 
formulation  ensures  a  highly  parallel  scheme  as  all  communication  is  local  between 
elements  sharing  edges  in  two-dimensions  or  faces  in  three-dimensions.  Finally,  the 
locality  of  the  approximation  makes  it  straightforward  to  extend  the  scheme  to 
include  support  for  different  orders  of  approximations  in  different  elements,  non- 
conforming  or  different  types  of  elements,  e.g.,  triangles  and  quadrilaterals. 

Before  discussing  the  numerical  flux,  let  us  note  that  a  mathematically  equivalent 
but  numerically  different  formulation  of  Eq.(5)  can  be  obtained  by  an  integration 
by  parts  to  recover 


-  Fn  •  V  -  SN^j  Lk(x )  dx  =  -j  Lk(x)F*N  dx  . 


(6) 


This  can  be  recognized  as  the  classical  discontinuous  Galerkin  method  for  conser¬ 
vation  laws  [1,  2,  4,  18].  To  distinguish  between  the  two  formulations,  we  shall 
refer  to  Eq.(5)  as  the  divergence  form  and  the  more  familiar  one,  Eq.(6),  as  Green’s 
form.  Other  terms  often  used  to  describe  these  formulations  are  the  strong  and  the 
weak  form,  respectively. 

The  numerical  flux,  F*N,  is  the  part  of  the  formulation  that  allows  information 
to  be  passed  between  the  individual  elements,  the  union  of  which  forms  S.  The 
discontinuous  formulation  implies  that  the  solution  at  an  interface  is  non-unique 
and  we  must  ensure  that  a  unique  solution  be  identified  and  passed  to  both  elements 
in  a  way  consistent  with  the  dynamics  of  the  problem. 

This  is  a  situation  similar  to  a  classical  finite  volume  formulation  to  which  Eq.(6) 
reduces  for  the  lowest  order  elements.  Thus,  we  can  borrow  from  the  extensive 
literature  devoted  to  the  development  and  analysis  of  numerical  fluxes  within  the 
context  of  finite  volume  methods,  e.g.,  upwinding  by  linearization  or  approximate 
Riemann  solvers  such  as  Roe  [23],  Engquist-Osher  [6]  or  Van  Leer  [28]  fluxes. 

For  simplicity  and  generality  we  use  subsequently  the  simple  Lax-Friedrichs  flux 
of  the  form 


*  FN(qN)  +  FN(pN)  |A| 
fn  = - ^ - Y  (Pv  —  Qn)  » 


where  qN  refers  to  the  local  computed  solution  and  pN  refers  to  the  solution  in 
the  neighboring  element.  The  dissipative  term  is  scaled  by  |A|  which  represents  the 
maximum  (local)  eigenvalue  of  the  flux  Jacobian 


0 
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where 


U  =  h  ■  u  . 

The  eigenvalues  of  the  flux  Jacobian  are  A  =  [U,U,U  +  y/ip,  U  —  ^JTp\ T ,  such  that 

|A|  =  |^|  +  v^  , 

is  the  maximum  wave  speed  of  the  shallow  water  equations  entering  the  Lax- 
Friedrichs  flux.  While  it  is  well  known  that  the  use  of  a  Lax-Friedrichs  flux  in  a 
classical  finite  volume  formulation  leads  to  a  very  dissipative  scheme,  this  is  much 
less  of  a  problem  in  a  high-order  formulation  where  the  quality  of  the  numerical 
flux  is  less  critical  [5]. 

To  simplify  matters  further  for  both  formulations,  let  us  introduce 


=  /  Li{x)Lk(x)dx  ,  T>ik=  /  Li(x)V Lk(x)  dx  ,  (7) 

JD  JD 

as  the  mass  matrix  and  the  differentiation  operator,  respectively.  Note  that  D  = 
[D^D^D1]  is  a  vector  of  matrices  corresponding  to  a  discrete  gradient  oeprator. 
To  account  for  the  source  we  define 


Sife  =  ^ 


(JV^)2 


R 2  ^ 

m= 1 


Xr. 


Lm  0 x)Li{x)Lk{x )  dx 


M&  = 


Li{x)Lk{x)S7ips{x)  dx 


As  D,  both  S  and  Mv  are  3- vectors  of  matrices. 

Finally  we  introduce  the  operator  associated  with  the  surface  integral  as 


F/fc  (p  L[ (^x^Lk (iu)  dx  , 

JSD 

where  l  includes  the  trace  of  nodes  on  the  face  of  D  only.  Denoting  the  local 
element-wise  grid  vector  of  the  geopotential  as  tpN  and,  correspondingly,  the  grid 
vector  for  the  three  momentum  components  as  <p>uN  and  qN  =  [ipN,(puN]T,  we 
can  now  express  the  semi-discrete  approximation  of  the  divergence  form,  Eq.(5),  as 


(I4  0  M)—qN  +  D  ■  F1  jv(<Zjv)  = 


0  0 

Sx 


Qn 


(8) 


0  0 
0  /rl3  <g>  M 


0 

xN  _ 

+ 14  <S>  F  [.Fjv(<7jv) 


where  Ir  is  a  rank-r  identity  matrix  and  the  Lagrange  multipliers  are  contained 
in  the  diagonal  matrix,  p  =  diag[/ii, ...,  //(n+i)2].  Finally  we  have  introduced  the 
arrays  of  nodal  physical  coordinates,  xjy  =  [xjv,  Un,  zn]T ■ 

Similarly,  we  can  express  the  Green’s  form,  Eq.(6),  using  the  above  notation  as 


(I4  <8  M)^-<Lv  -  •  FN(qN)  = 


0  0 

Sx 


Qn 


(9) 


0  0 
0  11I3  8  M 


0 

xN  _ 

(I4®F)F* 
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In  both  cases  we  can  compute  the  Lagrange  multipliers  in  a  way  similar  to  the 
continuous  case,  Eq.(3),  i.e.,  by  requiring  that  the  discrete  momentum  is  pointwise 
normal  to  the  position  vector.  Note  that  the  effect  of  the  Coriolis  force  needs  to  be 
included  in  the  discrete  form  to  ensure  that  fluid  particles  remain  on  the  spherical 
shell. 

It  is  worth  noticing  that  both  schemes  given  above  are  fully  explicit  in  time, 
i.e.,  no  global  assembly  is  required  in  contrast  to  a  classical  finite  element/spectral 
element  scheme,  and  is  thus  parallel  by  construction.  The  only  overhead  associated 
with  the  above  formulation  is  the  additional  memory  required  to  store  the  multiple 
solutions  at  the  overlapping  regions  of  the  elements,  i.e.,  the  edges  (sides).  For 
high-order  elements  this  memory  overhead  is  clearly  negligible  as  the  number  of 
edge  to  volume  nodes  decreases  rapidly  as  the  order  increases. 


3.2.  Temporal  Integration  and  Stability 

The  explicit  nature  of  the  semi-discrete  formulation 


dqN 

dt 


B(<7Ar)  , 


where  B  signifies  the  operators  given  in  Eqs.(8)-(9),  makes  it  natural  to  use  a 
standard  explicit  Runge-Kutta  scheme.  This  yields 


Vfc  =  1, ...,  3  :  qkN+1  =  qkN  +  AtakB(qk)  , 
with  qkjfl  =  q^f  and 


nn+1  —  nn 

Qn  ~  Qn 


(qk)  , 


k= 1 


where 

Oil  =a2  =  ^,  0:3  =  1,  and  (h  =  34  =  L  P2  =  /?3  =  2  . 

The  time  step,  At,  is  chosen  in  order  to  ensure  a  stable  scheme  and  will  generally 
scale  like 


At  <  CFL  x  min  [|n  •  *|  +  ip^/x  •  *J 
xgD 

where  the  local  grid  distortion  is  measured  by 

(M  +  M  M  +  M  JM  +  j uA\ 

X  Ai)’Ar  Ai)’A(  Ai)j  ' 

Here  (A^,  Aij)  reflects  the  local  average  grid  size  and  u  =  (u,  v,  w)  the  local  velocity. 

Even  with  a  suitably  chosen  value  of  the  time-step  it  is  well  known  that  high- 
order  methods  are  prone  to  instabilities  due  to  the  nonlinear  mixing  and  lack  of 
dissipation,  see  e.g.  [10].  This  is  particularly  true  for  problems  with  marginally 
resolved  phenomena  where  the  nonlinear  mixing  of  the  solution  with  the  Gibbs 
oscillations  can  drive  the  instability. 
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The  standard  approach  to  avoid  this  instability  in  a  controlled  manner  is  through 
the  use  of  a  weak  high-order  filter  which  modifies  the  high  frequency  modes  without 
altering  the  well  resolved  low  frequency  modes.  As  has  been  shown  over  the  last 
decade  such  filtering  can  be  applied  without  sacrificing  spectral  accuracy  [10,  26]. 

We  shall  focus  on  the  filters  developed  by  Boyd  [3]  and  Vandeven  [29].  While 
these  filters  perform  well,  their  use  are  known  to  pose  difficulties  in  a  classical 
spectral  element  scheme  where  the  solution  is  required  to  be  continuous.  As  we 
shall  see  shortly,  however,  these  difficulties  vanish  in  the  current  formulation. 
Following  the  discussion  in  [3],  consider  the  state  variables,  qN,  as 

(N+ 1)2  N 

Qn=  <lN(€k)Lk(€)  =  , 

k—1  i,j= 0 

where  _£/*,(£)  is  the  Lagrange  polynomial  associated  with  the  points,  £fe,  Pj(£)  and 
Pj(rf)  the  Legendre  polynomials  in  £  and  r/,  respectively,  and  q ^  the  discrete  Leg¬ 
endre  expansion  coefficients  of  the  state  vector  computed  by  using  the  Legendre 
quadrature. 

The  filtering  approach  proposed  in  [3]  involves  the  weighted  sum 

N 

Qn  =  (!  -  v)Qn  +  aia^ijpi(0pj(v)  ,  (10) 

i,j= 0 

where  v  is  the  filter  weighting,  i.e.,  v  —  0  represents  no  filtering  and  v  =  1  full 
filtering.  Furthermore 


o%  = 


1  for  i  <  s 

^(jrh)  for  s<i<N 


is  the  filter  function  with  a  being  the  Boyd- Vandeven  filter  [3];  a3  is  defined  likewise. 
Note  that  the  filtering  is  performed  in  an  element  sense.  In  other  words,  each 
element  yields  its  very  own  set  of  filtered  values.  However,  if  using  a  spectral  element 
formulation  we  would  need  to  require  that  the  state  vector  be  continuous  across 
elements.  As  the  filtering  gives  different  values  for  grid  points  shared  by  elements 
some  form  of  weighting  is  typically  required  to  recover  continuity.  Since  we  advocate 
a  discontinuous  element  formulation,  however,  this  correction  is  unnecessary,  hence 
greatly  simplifying  the  inclusion  of  filters  into  the  schemes. 


4.  GENERATION  OF  ICOSAHEDRAL  GRIDS  ON  THE  SPHERE 

Contrary  to  the  more  traditional  solution  techniques  exploiting  spherical  har¬ 
monics  [11,  14,  16,  24,  25],  the  use  of  a  multi-element  formulation  allows  for  the 
use  of  any  type  of  grid,  e.g.,  not  only  unstructured  grids  but  non-conforming  grids 
as  well. 

The  generation  of  the  grid  on  the  sphere  is  a  challenging  problem  and  in  this 
section  we  describe  the  procedure  for  generation  of  general  high  order  icosahedral 
grid  proposed  in  [8,  9].  This  grid  is  derived  from  the  icosahedron  comprised  of  20 
equilateral  triangular  elements  and  12  grid  points. 
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To  construct  icosahedral  grids  we  consider  the  initial  icosahedron  and  subdivide 
each  of  the  initial  triangles  by  a  triangular  Lagrange  polynomial  of  order  nj.  Prior 
to  mapping  these  elements  onto  the  sphere  it  is  convenient  to  map  the  triangles 
onto  a  gnomonic  space.  The  most  unbiased  mapping  is  obtained  by  mapping  about 
the  centroid  of  the  triangles. 

Let  (Ac,  9C)  be  the  centroid  of  the  triangle  we  wish  to  map.  The  gnomonic 
mapping  is  then  given  by 

acos#sin(A  —  Ac) 

x  - ,  fll) 

sin  9C  sin  9  +  cos  9C  cos  9  cos(A  —  Ac)  ’ 

a  [cos  9C  sin  9  —  sin  9C  cos  9  cos(A  —  Ac)] 

^  sin  #c  sin#  +  cos  #c  cos  #cos(A  —  Ac) 

To  simplify  matters  a  bit,  we  first  apply  a  rotation  whereby  Eq.(ll)  becomes 


x  =  a  tan  A',  y  =  a  tan  9'  sec  A'  ,  (12) 

in  the  new  coordinate  system  with  the  coordinates  (A ,9)  located  at  (0,0).  The 
rotation  mapping  is  given  as 


cos#sm(A  —  Ac) 

A  =  arctan  - — - — ; — - - - - - 7- - — T 

sin  9C  sm  9  +  cos  9C  cos  9  cos(A  —  Ac) 

9'  =  arcsin  [cos  9C  sin  9  —  sin  9C  cos  9  cos(A  —  Ac)] . 


(13) 


This  approach  enables  the  construction  of  a  general  icosahedral  grid  defined  by 


iVj  =  10(n/-l)2  +  20(n/-l)  +  12  , 

Nj  =  2(Np  —  2)  ,  (14) 

Nj  =  3(Np  —  2)  , 

where  Nj ,Nj,  and  Nj  denote  the  number  of  points,  elements,  and  sides  com¬ 
prising  the  triangular  grid. 

Once  the  triangular  icosahedral  grid  is  constructed,  we  subdivide  each  triangular 
element  into  3  quadrilateral  elements.  Upon  dividing  the  triangles  into  quadri¬ 
laterals  one  can  construct  the  higher  order  collocation  points  inside  each  element 
resulting  in  a  quadrilateral  grid  with  the  following  properties 

Np  =  6(Np  —  2)N2  +  2  , 

Ne  =  6  (JVpr  -  2)  ,  (15) 

Ns  =  12{Np  —  2)  , 


where  Np ,  Ne,  and  Ns  denote  the  number  of  points,  elements  and  sides  comprising 
the  quadrilateral  grid,  and  N  is  the  polynomial  order  used  in  the  semi-discrete 
discontinuous  Galerkin  scheme  discussed  in  Sec.  3.1.2. 
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TABLE  I 


The  number  of  grid  points, 
grid  as  a 

elements,  and  sides 
function  of  nj  and  N . 

for  the 

icosahedral 

ni 

N 

nj  •  IV 

Np 

Ne 

N3 

1 

4 

4 

962 

60 

1 

6 

6 

2162 

60 

120 

1 

8 

8 

3842 

60 

120 

1 

12 

12 

8642 

60 

120 

1 

16 

16 

15362 

60 

120 

1 

24 

24 

34562 

60 

120 

1 

32 

32 

61442 

60 

120 

1 

64 

64 

245762 

60 

120 

4 

1 

4 

962 

960 

1920 

8 

1 

8 

3842 

3840 

7680 

16 

1 

16 

15362 

15360 

30720 

32 

1 

32 

61442 

61440 

122880 

64 

1 

64 

245762 

245760 

491520 

2 

2 

4 

962 

240 

4 

2 

8 

3842 

960 

1920 

8 

2 

16 

15362 

3840 

7680 

16 

2 

32 

61442 

15360 

30720 

32 

2 

64 

245762 

61440 

122880 

1 

4 

4 

962 

60 

m 

2 

4 

8 

3842 

240 

19 

4 

4 

16 

15362 

960 

1920 

8 

4 

32 

61442 

3840 

7680 

16 

4 

64 

245762 

15360 

30720 

Substituting  the  values  in  Eq.(14)  into  Eq.(15)  yields 


Np 

=  60(n/)2lV2  +  2  , 

(16) 

Ne 

=  60(n/)2  , 

(17) 

Ns 

=  120(n/)2  . 

(18) 

Table  I  provides  examples  of  grids  for  various  values  of  n/  and  N.  Examples  of 
corresponding  grids  for  m  —  1  and  N  =  4,  8, 16,  and  32  are  illustrated  in  Fig.  2. 

5.  RESULTS 

In  the  following  we  evaluate  the  performance  of  the  scheme  discussed  in  the 
previous  sections.  As  a  measure  of  the  error  we  use  the  normalized  L2  error 


I|9jvIIl2 


'JpPZcxact  ~Qn)  dx 
J | )  dexaf:t 
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FIG.  2.  An  icosahedral  grid  for  nj  =  1  and  a)  N  =  4,  b)  N  =  8,  c)  N  =  16,  and  d)  N  =  32. 


throughout  this  section.  Here  qN  represents  the  computed  conservation  variables 
and  qexact  the  exact  when  available.  The  global  error  is  computed  as  a  broken  norm 
using  the  local  quadratures. 

Six  test  cases  are  used  in  order  to  test  the  algorithms  and  form  a  framework  for 
comparison  between  the  two  formulations.  Cases  1,  2,  3,  5  and  6  correspond  to  the 
test  cases  given  in  [30] .  Case  4  has  been  used  as  a  test  case  for  the  shallow  water 
equations  in  [8,  9,  20,  21].  Cases  1,  2,  and  3  have  analytic  solutions  and  are  used  to 
evaluate  the  accuracy  of  the  discontinuous  Galerkin  method  quantitatively.  Cases 
4,  5,  and  6,  on  the  other  hand,  do  not  have  analytic  solutions  and  are  thus  used  to 
obtain  a  qualitative  assessment  of  the  accuracy  of  the  scheme. 

•  Case  1:  Steady-State  Advection.  This  case  concerns  the  solid  body  rotation  of 
a  cosine  wave.  It  only  tests  the  mass  equation  as  the  velocity  field  is  assumed  to 
remain  unchanged  throughout  the  computation.  The  cosine  wave  completes  one 
full  revolution  after  12  days  of  integration. 
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•  Case  2:  Global  Steady-State  Nonlinear  Zonal  Geostrophic  Flow.  This  case  is  a 
steady-state  solution  to  the  nonlinear  shallow  water  equations.  The  equations  are 
geostrophically  balanced  and  remain  so  for  the  duration  of  the  integration.  The 
velocity  field  thus  remains  constant  throughout  the  computation.  The  geopotential 
height  ip  undergoes  a  solid  body  rotation  but  since  the  initial  height  field  is  given 
as  a  constant  band,  the  solution  remains  the  same  throughout  the  time  integration. 
The  velocity  field  is  the  same  as  that  used  in  Case  1.  Williamson  et  al.  [30] 
recommend  that  the  error  be  computed  after  5  days  of  integration. 

•  Case  3:  Steady-State  Nonlinear  Zonal  Geostrophic  Flow  with  Compact  Support.. 
This  case  is  similar  to  Case  2  except  that  the  velocity  is  zero  everywhere  except  in 
a  very  small  isolated  region.  This  isolated  region,  or  jet,  encapsulates  the  flow  and 
limits  the  geopotential  height  field  to  remain  within  a  very  confined  circular  region. 
As  in  Case  2,  the  errors  are  computed  after  5  days. 

•  Case  Dancing  High-Low  Waves.  This  case  comes  from  [20]  and  is  not  an 
analytic  solution  to  the  shallow  water  equations.  The  initial  geopotential  height  is 
comprised  of  two  large  waves  with  the  left  wave  being  the  low  wave  and  the  right 
wave  being  the  high  wave,  when  viewed  from  the  north  pole.  The  waves  rotate 
clockwise  in  a  swirling  dance-like  fashion  so  that  after  10  days  of  integration,  the 
low  wave  is  again  on  the  left  and  the  high  wave  is  on  the  right. 

•  Case  5:  Zonal  Flow  over  an  Isolated  Mountain.  This  case  is  similar  to  Case  2 
except  that  a  mountain  has  been  included  on  the  sphere.  This  is  the  only  problem 
in  the  test  cases  studied  here  which  includes  topography.  The  mountain  is  conical 
in  shape  and  forces  the  zonal  flow  to  deflect  off  the  mountain.  Due  to  the  zonal 
flow  impinging  on  the  mountain,  wave  structures  form  and  propagate  around  the 
sphere.  Results  are  reported  for  a  10  day  integration  period. 

•  Case  6:  Rossby-Haurwitz  Wave.  Although  Rossby-Haurwitz  waves  are  not 
analytic  solutions  to  the  shallow  water  equations,  they  have  become  a  standard 
test  case.  In  a  non-divergent  barotropic  model,  the  initial  geopotential  height  field 
undergoes  a  solid  body  rotation  in  a  counterclockwise  direction  (when  viewed  from 
the  north  pole).  Results  are  reported  for  a  14  day  integration. 

In  the  following  we  shall  use  these  6  benchmarks  as  the  stick  against  which  to 
measure  and  compare  the  schemes  and  their  numerical  properties  such  as  robustness 
and  accuracy.  For  the  former  we  shall  discuss  the  impact  of  various  simplifications 
and  approximations  introduced  into  the  two  schemes,  Eqs.(8)-(9),  while  the  latter  is 
evaluated  by  comparison  with  exact  solutions  as  well  as  studies  of  a  more  qualitative 
nature. 


5.1.  Basic  Convergence  Tests 

Prior  to  that,  however,  it  seems  only  natural  to  discuss  and  illustrate  the  advan¬ 
tages  of  using  a  high-order  scheme  for  solving  the  shallow  water  equations.  For  this 
purpose  we  shall  consider  the  solution  of  Cases  1,  2,  and  3  using  the  divergence 
form  of  the  scheme,  Eq.(8),  only. 

Figure  3  shows  the  computed  mass  error  for  the  3  cases  using  the  following 
orders  of  accuracy:  N  =  1  (dashed),  N  =  2  (dotted),  N  =  4  (dashed-dotted),  and 
n/  =  1  (solid).  We  plot  the  mass  error  norm  as  a  function  of  the  product  n/N.  The 
N  =  1  results  are  obtained  by  using  linear  elements  but  increasing  n/  for  increasing 
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FIG.  3.  The  mass  error  norm  for  the  divergence  form  of  the  DGM  comparing  low  order 
schemes  ( N  <  4)  with  the  Nth  order  scheme  ( n /  =  1)  for  a)  Case  1,  b)  Case  2,  and  c)  Case  3  for 
a  one  day  integration. 


niN,  i.e. ,  it  corresponds  to  a  classic  element  refinement  known  as  /i-refinement.  In 
contrast,  the  nj  =  1  results  shown  in  the  plots  represent  our  JVth  order  scheme 
where  we  keep  the  number  of  elements  constant  (given  by  n/  =  1)  and  increase  N 
as  is  done  in  classic  high-order/spectral  convergence. 

The  results  in  Fig.  3  confirm  the  spectral  accuracy  of  the  scheme  for  all  three 
test  cases.  These  results  also  illustrate  well  the  advantages  offered  by  high-order 
schemes  over  low-order  schemes  in  terms  of  accuracy. 
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However,  it  is  crucial  to  understand  the  cost  incurred  by  this  additional  accuracy, 
e.g.,  if  the  high-order  scheme  is  prohibitively  expensive  for  levels  of  accuracy  of 
relevance  then  the  scheme  loses  much  of  its  appeal. 

In  order  to  measure  cost  versus  accuracy,  we  shall  consider  the  number  of  op¬ 
erations  associated  with  the  different  schemes.  The  simplest  reasonable  operation 
count  is  0(N5Ne  /2),  where  N  is  the  local  approximation  order  on  the  Ne  elements. 
This  estimate  is  obtained  by  including  the  Ne  evaluations  of  the  derivatives,  being 
0(N3Ne),  and  including  an  additional  0(N2\fWl)  work  to  account  for  the  number 
of  time-steps  required  to  advance  the  high-order  method.  Taking  the  expression  for 
Ne  from  Eq.(18),  this  yields  an  estimate  of  the  work  as  0(N5n3). 

The  main  interest  is  naturally  to  limit  the  number  of  operations  needed  to  achieve 
a  result  with  a  given  accuracy,  which,  on  the  other  hand  is  a  problem  dependent 
entity.  Let  us  first  consider  the  results  of  Case  1,  illustrated  in  Fig.  3a.  Taking 
an  accuracy  of  10-3  as  the  goal,  it  is  easily  found  from  Fig.  3a  that  the  4th  order 
scheme  is  the  one  with  the  least  operations,  about  half  that  of  the  2nd  order  scheme 
and  almost  5  times  less  than  the  first  order  scheme.  In  a  similar  fashion  one  finds 
for  Case  2  and  an  accuracy  goal  of  10-6  that  the  Nth  order  scheme  seems  most 
efficient  while  the  1st  and  2nd  order  methods  are  prohibitively  expensive.  For  Case 
3  and  an  accuracy  goal  of  ICC3  we  again  find  that  the  Nth  order  scheme  appears 
to  be  most  efficient. 

These  results  provide  only  guidelines  but  they  do  confirm  the  advantages  in 
using  a  high-order  scheme  over  lower  accuracy  methods.  Whether  one  should  use 
the  highest  order  approximation  possible  or  rather  limit  the  order  and  refine  the 
element  size  may  well  be  problem  dependent.  However,  as  has  also  been  found  in 
other  similar  studies  [9,  12,  13,  15,  26]  it  is  generally  advantageous  to  use  a  moderate 
order  of  approximation,  N  =  8  —  16,  and  refine  the  element  size  accordingly  to 
achieve  a  practical  level  of  accuracy.  Only  for  problems  requiring  very  high  accuracy 
or  long  time  integration  can  one  benefit  from  using  a  very  high  order  scheme,  i.e., 
N  >  16.  The  results  obtained  here  support  the  validity  of  these  guidelines. 


5.2.  The  Divergence  Form  versus  the  Green’s  Form 

As  a  second  test  we  shall  evaluate  the  differences  between  the  divergence  form, 
Eq.(8),  and  the  Green’s  formulation,  Eq.(9),  subject  to  various  approximations 
and  simplifications.  As  a  basis  for  this  comparison  we  shall  again  use  Cases  1,  2, 
and  3,  with  the  comparisons  including  the  impact  of  using  full  or  lumped  diagonal 
mass  matrices,  and  the  use  of  exact  or  inexact  integration  to  compute  the  discrete 
operators.  Although  these  differences  appear  as  being  small  we  shall  see  that  the 
performance  of  the  divergence  and  Green’s  formulations  can  be  quite  different  as  a 
consequence  of  these  differences. 


5.2.1.  Full  versus  Diagonal  Mass  Matrix 

One  of  the  most  immediate  ways  of  improving  the  efficiency  of  the  discontinuous 
Galerkin  method  is  to  approximate  the  mass  matrix,  M,  Eq.(7),  by  a  diagonal  (or 
lumped)  form.  We  form  the  diagonal  approximation  of  the  mass  matrix  by  simply 
summing  all  of  the  entries  of  each  row  and  storing  the  sum  in  the  main  diagonal. 
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FIG.  4.  Evaluation  of  the  ip  error  induced  by  using  diagonal  (square)  and  full  (triangle) 
mass  matrices,  M,  in  the  DGM  formulations.  The  error  norms  are  illustrated  for  the  divergence 
form  (solid  line)  and  the  Green’s  form  (dashed  line)  for  a)  Case  1,  b)  Case  2,  and  c)  Case  3.  All 
results  correspond  to  one  day  integrations. 


Figure  4  shows  the  errors  obtained  using  the  full  and  diagonal  mass  matrices  for 
the  divergence  and  Green’s  forms,  Eqs.(8)-(9),  on  the  nj  =  1  grid  for  Cases  1,  2, 
and  3. 

The  results  illustrate  that  the  divergence  form  is  generally  more  robust  towards 
the  approximation  induced  by  the  diagonal  mass  matrix,  exemplified  most  clearly 
by  the  solution  of  Case  2.  Note,  however,  that  the  disparity  in  accuracy  between 
the  full  and  diagonal  forms  decreases  for  all  three  cases  as  the  order  of  the  basis 
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functions  is  increased.  The  level  off  of  the  error  decay  in  Case  2  we  attribute  to 
finite  precision  effects  in  combination  with  doing  the  computations  in  dimensional 
variables.  This  effect,  however,  does  not  impact  the  conclusions. 

5.2.2.  Exact  versus  Inexact  Integration 

The  evaluation  of  the  discrete  operators,  e.g.,  Eq.(7),  requires  the  computations 
of  integrals.  Despite  the  fact  that  we  have  the  solution  defined  on  Legenclre-Gauss- 
Lobatto  quadrature  points,  the  curvilinear  geometry,  reflected  by  the  transforma¬ 
tion  Jacobian,  and  the  terms  associated  with  the  Coriolis  force  will  create  polynomi¬ 
als  of  higher  order  and  over-integration  would  be  needed  to  integrated  these  terms 
exactly.  Indeed,  to  integrate  all  of  the  matrices  exactly  requires  Q  =  ( cN  +  l)/2 
quadrature  points,  where  c  is  an  integer  constant  denoting  the  factor  of  the  maxi¬ 
mum  order  matrix.  For  the  spherical  shallow  water  equations  c  =  4,  coming  from 
the  Coriolis  term.  For  highly  skewed  elements,  however,  this  factor  increases  to  6 
because  of  the  impact  of  the  transformation  Jacobians. 

The  question  to  address  here  is  what  is  the  effect  of  employing  inexact  in¬ 
tegration  only,  as  is  done  traditionally  in  continuous  spectral  element  schemes 
[7,  15,  19,  26].  Here  we  shall  simply  use  the  straightforward  quadrature  associ¬ 
ated  with  the  Legendre-Gauss-Lobatto  nodes,  to  evaluate  the  inner  products  and, 
thus,  the  discrete  operators. 

Figure  5  shows  the  mass  error  obtained  using  exact  and  inexact  integration  for 
the  matrices  of  the  divergence  and  the  Green’s  forms  on  the  nj  =  1  grid  for  Cases 
1,  2,  and  3.  Both  the  divergence  and  the  Green’s  forms  clearly  loose  accuracy 
when  using  inexact  integration.  Similar  to  the  discussion  for  the  lumped  mass 
matrix,  however,  we  observe  that  the  divergence  form  appears  to  be  more  robust 
towards  such  approximations,  i.e.,  the  results  obtained  using  the  divergence  form  is 
always  less  affected  by  the  approximations  than  the  Green’s  form.  Furthermore,  we 
observe  that  this  increased  error  level  for  the  Green’s  form  seems  to  persist  for  large 
orders  of  approximation  in  Case  2.  In  other  words,  the  accuracy  disparity  between 
the  exact  and  inexact  integration  of  the  Green’s  form  appears  to  not  decrease  for 
increased  basis  function  order  as  it  does  for  the  divergence  form. 

5.2.3.  Divergence  versus  Green’s  Form 

The  results  in  Figs.  4  and  5  indicate  that  the  Green’s  form,  Eq.(9),  is  superior  in 
accuracy  to  the  divergence  form,  Eq.(8),  of  the  discontinuous  Galerkin  formulation 
only  if  everything  is  done  exactly,  i.e.,  full  mass  matrices  and  over-integration  to 
evaluate  the  discrete  operators  exactly.  However,  for  all  cases,  the  divergence  form 
performs  almost  as  well  even  when  doing  all  operations  exactly  and  it  was  found  to 
be  more  robust  to  performance  enhancing  approximations  such  as  diagonal  mass 
matrices  and  inexact  integration. 

Based  on  these  results  we  shall  use  the  divergence  form  with  inexact  integration 
and  diagonal  mass  matrices  for  the  remainder  of  the  paper  as  these  approximations 
have  a  negligible  impact  on  the  accuracy.  It  should  be  noted,  however,  that  the 
problems  we  have  considered  and  on  which  we  have  based  this  choice  are  limited 
and  it  is  not  clear  which  formulation  to  choose  for  problems  of  a  more  general 
nature,  e.g.,  problems  with  discontinuous  solutions.  We  hope  to  address  this  in 
future  work. 


NODAL  DGM  FOR  THE  SPHERICAL  SHALLOW  WATER  EQUATIONS 


19 


FIG.  5.  Evaluation  of  the  ip  error  induced  by  using  exact  (triangle)  and  inexact  (square) 
integration  for  the  computation  of  the  operators  in  the  DGM  formulations.  The  error  norms  are 
illustrated  for  the  divergence  form  (solid  line)  and  the  Green’s  form  (dashed  line)  for  a)  Case  1, 
b)  Case  2,  and  c)  Case  3.  All  results  correspond  to  one  day  integrations. 


5.3.  Convergence  and  Stability  Study 

In  this  section  we  show  convergence  results  for  all  six  cases.  Cases  1,  2,  and 
3  have  analytic  solutions  and  so  we  use  these  cases  to  judge  the  accuracy  of  our 
results  quantitatively.  Cases  4,  5,  and  6,  on  the  other  hand,  are  only  used  to  judge 
the  convergence  of  the  scheme  qualitatively,  as  we  do  not  have  analytic  solutions 
to  these  cases.  All  six  test  cases  are  integrated  for  long  periods  of  time  in  order  to 
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confirm  the  stability  of  the  scheme.  To  ensure  stability  we  apply  the  Boyd-Vandeven 
filter  after  every  10  time  steps. 

5.3.1.  Quantitative  Study 

Figure  6  shows  the  mass  error  norm  as  a  function  of  the  approximation  order, 
N,  for  Cases  1,  2,  and  3.  The  results  for  Case  1  are  for  a  12  day  integration 
which  corresponds  to  a  complete  revolution  of  the  cosine  wave  around  the  sphere. 
For  Cases  2  and  3,  the  results  are  for  5  day  integrations  which  is  the  time  frame 
recommended  in  [30]. 

For  all  cases  we  observe  exponential  convergence  rates  until  the  effect  of  the 
finite  precision  impacts  the  solution  accuracy.  To  confirm  that  the  leveling  off  is 
due  to  finite  precision  effects  we  have  performed  studies  at  a  decreased  time-step, 
confirming  that  the  dominating  error  source  is  spatial  and  not  temporal. 


FIG.  6.  The  ip  error  as  a  function  of  the  spatial  approximation  order,  N ,  for  Case  1 
(dashed),  Case  2  (dotted),  and  Case  3  (solid)  after  12,  5,  and  5  days  of  integrations,  respectively. 
The  72/  —  1  grid  is  used. 


5.3.2.  Qualitative  Study 

In  the  previous  section  we  tested  cases  with  analytic  solutions  against  which  to 
compare  and  judge  the  accuracy  of  the  schemes.  However,  the  more  relevant  Cases 
4,  5,  and  6,  do  not  allow  such  simple  analytic  solutions.  Instead  we  shall  evaluate 
the  convergence  characteristics  of  the  scheme  qualitatively  by  running  the  cases 
using  different  resolutions. 

In  Figs.  7,  8,  and  9  the  left  panels  show  the  contours  of  the  mass,  the  middle 
panels  show  the  zonal  velocity  component,  and  the  right  panels  show  the  meridional 
velocity  component.  The  zonal  velocity  component  is  the  u  component  in  spherical 
coordinates  which  is  associated  with  the  longitude  A  and  is  positive  when  traveling 
west  to  east  in  a  direction  parallel  to  the  Equator.  The  meridional  velocity  com¬ 
ponent  is  the  v  component  in  spherical  coordinates  which  is  associated  with  the 
latitude  6  and  is  positive  when  traveling  from  the  South  Pole  to  the  North  Pole. 

Figure  7  shows  the  results  for  Case  4  using  the  ni  =  1  grid  with  various  orders 
of  approximation,  N ,  for  a  10  day  integration  and  confirm  that  the  wave  structures 
remain  the  same  for  all  three  values  of  N.  However  the  wave  pattern  clearly  becomes 
better  resolved  when  increasing  N .  This  is  most  obvious  near  the  North  Pole  (center 
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a) 


b) 


c) 


FIG.  7.  Case  4.  Contours  of  the  mass  (left),  u- velocity  (middle),  and  v- velocity  (right)  on 
grid  nj  =  1  and  a)  N  =  8,  b)  N  =  16,  and  c)  N  =  32  after  a  10  day  integration  viewed  from 
(A,  6)  =  (0,  90). 


of  the  plots)  for  the  velocity  components.  The  contours  are  a  bit  jagged  for  N  =  8 
but  become  smoother  for  N  =  16.  Finally,  for  N  =  32  we  see  that  the  same 
patterns  exist  as  in  the  N  =  8  and  TV  =  16  contours  but  the  curves  are  much 
smoother  throughout,  indicating  a  converged  result. 

Figure  8  shows  the  results  of  Case  5  for  a  10  day  integration  using  N  =  8, 16, 
and  32,  respectively.  The  contours  are  shown  from  the  viewpoint  (A,0)  =  (180,0), 
where  the  peak  of  the  mountain  resides  at  (A,  6)  =  (180,  30).  The  flow  impinging  on 
the  mountain  causes  the  resulting  wave  structures  which  we  see  in  the  figure  around 
the  mid-latitudes  (9  =  ±45).  The  contours  for  N  =  8  are  again  a  little  jagged  - 
not  just  for  the  mass  contours  but  for  the  velocity  components  as  well.  Increasing 
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FIG.  8.  Case  5.  Contours  of  the  mass  (left),  u-velocity  (middle),  and  v-velocity  (right)  on 
grid  nj  =  1  and  a)  IV  =  3,  b)  A  =  16,  and  c)  N  =  32  after  a  10  day  integration  viewed  from 
(A,  6)  =  (180,0). 


N  to  16  removes  much  of  the  jaggedness  visible  in  the  contours  and  increasing  N 
further  to  32  smoothen  the  contours  completely  resulting  in  clearly  visible  cohesive 
wave  structures. 

Figure  9  shows  the  results  of  Case  6  after  a  14  day  integration  at  different  resolu¬ 
tions  with  the  contours  being  shown  from  a  viewpoint  corresponding  to  the  North 
Pole  (A ,0)  =  (0,90).  In  contrast  to  the  previous  two  cases,  the  results  for  N  —  8 
show  that  the  wave  structures  are  beginning  to  break  down  at  this  resolution  which 
is  insufficient  to  support  the  dynamics  of  the  system.  This  breakdown  of  the  wave 
structures  is  most  obvious  by  looking  at  the  mass  and  velocity  field.  Increasing  N 
to  16  results  in  a  dramatic  improvement  with  the  wave  structures  becoming  more 
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cohesive.  This  is  particularly  noticeable  in  the  mass  contours  where  we  see  that  the 
rotated  Greek  cross  pattern  remains  intact  instead  of  breaking  off  into  distinct  blobs 
as  in  Fig.  9a.  Increasing  TV  further  to  32  significantly  improves  the  wave  patterns. 
The  mass  contours  for  the  TV  =  16  scheme  are  beginning  to  reveal  a  semblance  of 
break-off  of  the  wave  pattern.  This  can  be  seen  by  looking  at  the  bottom  left  tip 
of  the  cross.  For  TV  =  32  the  wave  pattern  remains  completely  intact  without  any 
indication  of  break-off. 


FIG.  9.  Case  6.  Contours  of  the  mass  (left),  u-velocity  (middle),  and  v-velocity  (right)  on 
grid  iif  =  1  and  a)  N  =  8,  b)  N  =  16,  and  c)  N  =  32  after  a  14  day  integration  viewed  from 
(A,  9)  =  (0,  90). 


The  results  for  these  three  cases  show  that  increasing  the  order  of  the  approxi¬ 
mation,  TV,  either  improves  the  smoothness  of  the  contour  curves  or  allows  for  the 
resolution  of  finer  scale  waves.  This  is  particularly  noticeable  in  Cases  5  and  6.  In 
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Case  5,  the  velocity  components  exhibit  very  localized  wave  formations  that  are 
extremely  well  resolved  by  our  nodal  IVth  order  DGM  scheme.  Case  6  illustrates 
the  breakdown  of  the  wave  structure  if  insufficient  grid  resolution  is  used. 

6.  CONCLUSIONS 

The  objective  of  this  paper  has  been  to  present  the  formulation  and  verification  of 
a  high-order  accurate  nodal  discontinuous  Galerkin  formulation  for  the  solution  of 
the  spherical  shallow  water  equations.  Curvilinear  quadrilateral  elements  are  used 
to  cover  the  sphere,  using  an  icosahedral  grid  as  the  basis  for  the  grid  generation, 
and  the  equations  are  solved  in  Cartesian  form  to  avoid  problems  with  coordinate 
singularities.  On  each  curvilinear  element  the  solutions  are  represented  by  Lagrange 
polynomials  in  a  purely  nodal  form,  i.e.,  the  solutions  are  given  on  quadrature 
points  and  operations  such  as  differentiation  and  integration  become  matrix-matrix 
operations.  The  equations  are  satisfied  in  a  discontinuous  element  form  with  the 
element  continuity  being  imposed  only  weakly.  This  decouples  all  elements  and 
makes  the  formulation  highly  parallel  as  well  as  well-suited  for  adaptive  solution 
techniques  as  no  constraints  on  element  conformity  is  needed. 

The  accuracy  of  the  scheme,  given  in  two  mathematically  equivalent  but  com¬ 
putationally  different  forms,  has  been  illustrated  by  considering  the  solution  of  the 
standard  set  of  benchmarks  proposed  in  [30].  The  results  confirm  the  expected 
high-order/spectral  accuracy  and  illustrate  the  advantages  of  using  such  methods 
to  efficiently  solve  geophysical  flow  problems. 

To  understand,  however,  whether  the  approach  proposed  here  provides  a  realistic 
alternative  to  existing  methods  based  on  spherical  harmonics  or  continuous  spectral 
element  methods  requires  extensive  further  validation.  The  inherent  properties 
of  the  proposed  technique,  e.g.,  high  parallel  efficiency,  high-order  accuracy  and 
support  for  adaptive,  non-conforming  solution  techniques,  are  sufficient  to  warrant 
such  exhaustive  studies  and  we  hope  to  report  on  such  results  in  the  near  future. 
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